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An  analytical  and  finite  element  model  of  a  single,  anode  supported  solid  oxide  fuel  cell  has  been  devel¬ 
oped  in  order  to  predict  the  stress  in  ceramic  components  subjected  to  an  idealised  operating  duty  cycle 
representing  cooling  from  sintering,  warming  to  a  uniform  temperature  of  800  °C  where  anode  chemical 
reduction  takes  place,  operation  at  low,  medium  and  high  power  and  finally  cooling  to  room  temperature. 

An  Abaqus™  finite  element  model  used  the  temperature  distribution  predicted  by  a  computational 
fluid  dynamics  model  at  low,  medium  and  high  power  to  solve  for  the  stress  distribution  throughout 
the  duty  cycle.  The  finite  element  model  included  the  effects  of  thermal  expansion,  residual  stress  from 
manufacture,  material  properties  changes  due  to  chemical  reduction  of  the  anode  and  visco-plastic  creep. 
The  level  of  stress  relaxation  predicted  by  the  finite  element  model  is  significant  at  SOFC  operating 
temperatures  and  timescales  of  several  thousand  hours. 

An  analytical  model  of  the  stress  distribution  in  thin  multilayer  plates  where  the  layers  have  different 
coefficients  of  thermal  expansion  was  developed  to  cross  check  the  finite  element  model.  In  the  analytical 
model  the  multilayer  plate  is  either  free  to  bend  or  constrained  to  remain  flat.  The  maximum  principal 
stresses  predicted  by  the  analytical  and  finite  element  models  were  found  to  agree  to  within  4%. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  typically  operate  at  a  temperature  in  excess 
of  800  °C  and  are  usually  stacked  and  clamped  together  to  pro¬ 
duce  a  usable  voltage.  During  operation  the  individual  cells  are 
subjected  to  high  stresses  caused  primarily  by  differential  thermal 
expansion,  mechanical  constraints  and  redox  reactions.  This  paper 
describes  the  development  of  a  thermal  stress  model  of  the  solid 
oxide  fuel  cell  that  takes  account  of  the  evolution  of  the  stress  state 
of  the  SOFC  from  sintering  (manufacture),  cooling  to  room  tem¬ 
perature,  uniform  heating  to  operating  temperature,  reduction  of 
the  NiO  anode  and  through  a  simplified  electricity  generating  duty 
cycle  before  final  cool-down,  representing  a  very  simplified  SOFC 
life-cycle.  The  temperature  distribution  across  the  fuel  cell  mem¬ 
brane  during  the  electricity-generating  phase  is  calculated  using  a 
detailed  CFD  based  electrochemical  model  developed  by  the  author 
in  previous  work  [5]. 

A  number  of  researchers  have  analysed  the  stress  state  of  SOFC 
layers  during  parts  of  the  SOFC  duty  cycle  such  as;  manufactur¬ 
ing  (residual  stress)  [24,3,7];  differential  thermal  expansion  [14]; 
temperature  gradients  [21,25];  volume  changes  on  reduction  or 
oxidation  [16];  external  mechanical  loading  or  restraint;  and  oxy¬ 
gen  activity  gradients  [2].  Chokshi  [4]  and  Morales-Rodriguez  et  al. 
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[15]  investigated  visco-plastic  creep  in  SOFC  materials,  which  acts 
to  relieve  stresses  over  long  timescales  compared  to  the  thermo¬ 
elastic  mechanisms.  There  is  also  some  evidence  [19]  that  allowing 
redox  cycling  to  occur,  where  the  anode  is  allowed  to  oxidise  on 
cooling  and  subsequently  must  be  reduced  again  when  the  cell 
is  heated  to  generate  electricity,  may  induce  additional  stresses. 
Although  a  single  anode  reduction  event  is  modelled  in  this  work 
multiple  redox  cycles  are  not  analysed. 

Many  of  the  stress  analyses  currently  presented  in  the  litera¬ 
ture  are  not  validated,  indeed  it  is  very  difficult  to  do  so.  Therefore 
a  general  analytical  solution  is  presented  for  the  stresses  due  to 
differential  thermal  expansion  in  a  thin  multi-layer  plate  either 
constrained  to  remain  flat  or  free  to  bend.  In  this  work  the  finite 
element  model  is  compared  to  the  analytical  solution  to  provide 
some  degree  of  validation. 

Ceramics  are  brittle  materials  and  as  such  their  bulk  failure 
strength  is  dependent  on  the  size  and  orientation  of  flaws  and  inclu¬ 
sions  that  act  as  stress  raisers.  This  has  two  effects:  firstly  it  gives 
a  large  spread  of  failure  strength  values;  secondly  it  gives  the  fail¬ 
ure  strength  a  volume  dependence  as  a  larger  piece  of  material  can 
contain  a  larger  stress  raiser  than  a  smaller  piece  [23].  Therefore  it 
is  more  usual  to  talk  in  terms  of  a  probability  of  failure  of  a  ceramic 
structure  rather  than  a  specific  stress  at  which  failure  occurs.  This 
makes  discussion  of  whether  a  certain  stress  level  is  critical  or  not 
for  a  ceramic  structure  difficult  without  an  in-depth  analysis  of  the 
probability  of  failure,  which  is  beyond  the  scope  of  this  paper.  A  first 
order  approximation  is  to  compare  the  stress  level  experienced  by 
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Top  Side  Bottom 


Fig.  1.  FEA  model  computational  domain. 


Nomenclature 

Q  activation  energy  (J  mol-1) 

h  enthalpy  or  cumulative  laminated  plate  thickness 

(J  mol-1  or  m) 

t  time  or  laminate  thickness  (s  or  m) 

c  uniform  strain  component  (-) 

tb  bending  axis  position  in  direction  x3  (m) 


the  ceramic  to  a  characteristic  strength,  defined  as  the  stress  at 
which  50%  of  samples  fail.  Caution  should  still  be  exercised  unless 
this  data  is  obtained  for  the  structure  in  question,  as  data  in  the 
literature  may  not  be  from  specimens  of  the  same  volume. 

2.2.  Model  geometry 

The  computational  domain  represents  a  typical  anode  sup¬ 
ported  SOFC  and  is  based  on  the  geometry  of  the  IndecASCl  cell 
shown  in  Fig.  1.  The  cell  dimensions  are  given  in  Table  1.  In  this 
work  a  single  fuel  cell  is  considered  to  be  bonded  to  a  test  housing 
over  the  area  of  surface  as  shown  in  Fig.  2.  The  test  housing  has 
gas  channels  cut  into  each  side  to  allow  air  and  hydrogen  to  flow 
over  the  cathode  and  anode  respectively. 

2.2.  Duty  cycle  definition 

A  duty  cycle  was  defined  to  represent  a  simplified  ‘real  world’ 
SOFC  operating  duty  cycle.  It  reflects  the  temperature  cycling 
and  mechanical  constraints  that  a  fuel  cell  membrane  might  be 


Table  1 

Dimensions  of  IndecASCl  cell. 


X\ <a  =  50  mm 

X2,a  =  10  mm 

hi  =  600  jxm 

Xib  =  60  mm 

x2,b  =  110  mm 

h2  =  610  |jim 

X\ iC  =  55  mm 

x2,c  =  120  mm 

h  3  =  616  jxm 

x2,d  =  5  mm 

h4  =  626  |xm 

x2,e  =  115  mm 

h5  =  656  |xm 

subjected  to  from  the  point  of  manufacture  onward.  Fig.  3  shows 
how  the  duty  cycle  is  defined  over  time. 

Initially  the  membrane  is  at  sintering  temperature,  then  it  is 
cooled  at  2°Cmin-1  to  room  temperature,  a  total  period  of  lOh. 
Normally  the  cell  would  then  be  stored  until  it  is  assembled  into  a 
stack  but  in  this  work  any  changes  in  stress  over  time  at  room  tem¬ 
perature  are  considered  negligible.  This  time  period  is  considered 
to  be  zero  in  the  analysis  but  still  indicated  on  the  duty  cycle  figure 
as  a  horizontal  line.  When  the  cell  gets  assembled  into  a  stack  or 
test  housing  the  mechanical  boundary  conditions  are  imposed,  in 
this  case  a  5  mm  wide  strip  around  the  outside  of  the  anode  layer 
is  restrained  in  the  Z-direction  to  simulate  bonding  the  cell  into 
the  housing.  The  cell  is  then  heated  in  a  furnace  to  800  °C  over  a 
period  of  nearly  7h  corresponding  to  a  rate  of  2°Cmin-1,  during 
this  time  the  temperature  distribution  in  the  cell  is  considered  to 
be  uniform  as  the  thin  cell  is  in  quasi-static  thermal  equilibrium. 
At  800  °C  hydrogen  is  gradually  introduced  to  the  fuel  stream  to 
chemically  reduce  the  NiO/YSZ  anode  to  Ni/YSZ  over  a  period  of 


Fig.  2.  Fluid  domain  of  SOFC  in  a  test  housing. 
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Table  2 

Material  properties  of  solid,  fully  dense  SOFC  materials. 


Temp 

Electrolyte 

YSZ 

Cathode 

LSM 

Anode 

NiO/YSZ 

Eo 

298  K 

190  [2] 

110  [8] 

207.2  [20] 

(GPa) 

1073  K 

157  [2] 

- 

171.2a 

1273  K 

- 

118  [8] 

- 

Vo 

298  K 

0.308  [2] 

0.36  [8] 

0.328  [20] 

1073  K 

0.313  [2] 

- 

0.333a 

1273  K 

- 

0.36  [8] 

- 

a 

298  K 

7.6  [11] 

9.8  [26] 

11.7  [26] 

(x  10-6) 

1273  K 

10.5  [11] 

11.8  [26] 

12.5  [26] 

a  Calculated  here  assuming  same  temperature  dependence  as  pure  YSZ. 


Time  (hrs) 


Fig.  3.  Simplified  SOFC  duty  cycle. 

3  h  to  produce  an  active  anode.  At  20  h  a  sequence  of  running  the 
fuel  cell  at  low  power  (2000  h),  high  power  (2000  h)  and  medium 
power  (2000  h)  is  implemented  before  cooling  the  fuel  cell  to  room 
temperature  at  a  rate  of  -2  °C  min-1  between  6020  and  6027  h. 

In  both  the  finite  element  and  analytical  models  duty  cycle 
the  temperature  boundary  conditions  are  interpolated  from  a  CFD 
model,  in  this  case  low,  medium  and  high  power  are  derived  from 
the  temperature  distribution  predicted  by  the  CFD  model  at  70, 1 60 
and  300  mA cm-2,  as  described  in  Section  4.6. 

2.  Material  properties 

In  this  work  the  materials  are  assumed  to  be  homogeneous 
and  isotropic.  The  following  sections  present  the  elastic  and  visco¬ 
plastic  material  properties. 

2.2.  Thermo-elastic  material  properties 

The  thermo-elastic  material  properties  available  in  the  litera¬ 
ture  for  solid,  fully  dense  material  samples  relevant  to  SOFCs  are 


summarised  in  Table  2  as  a  function  of  temperature.  Where  there 
is  a  variation  of  material  properties  such  as  a  and  E  with  temper¬ 
ature  a  linear  interpolation  is  assumed  between  the  data  points 
quoted  from  the  literature. 

The  elastic  modulus  and  Poisson’s  ratio  of  the  porous  cathode 
(LSM),  anode  (Ni/YSZ  or  NiO/YSZ)  and  a  porous  magnesia  spinel 
(MMA)  material  used  as  a  support  structure  in  the  Rolls  Royce  IP- 
SOFC,  have  all  been  shown  to  be  strongly  dependent  on  porosity 
[20,1 6,8].  This  is  particularly  relevant  for  the  material  properties  of 
the  anode,  where  the  volume  porosity  of  the  porous  anode  structure 
increases  significantly  on  reduction  from  NiO/YSZ  to  Ni/YSZ,  with  a 
corresponding  change  in  material  properties.  From  a  separate  study 
of  the  mechanical  properties  of  the  Ni-NiO/YSZ  cermet  as  a  function 
of  percentage  of  NiO  reduced  Radovic  and  Lara-Curzio  [16]  con¬ 
cluded  that  the  changes  in  bulk  material  properties  on  reduction 
were  due  predominantly  to  the  change  in  bulk  volume  porosity. 
There  have  been  several  models  proposed  to  account  for  this  poros¬ 
ity  dependence,  the  more  frequently  quoted  are;  the  exponential 
or  Minimum  Solid  Area  (MSA)  model  used  by  Spriggs  [22],  Knudsen 
[13]  and  Atkinson  [20]  based  on  the  reduction  in  load  bearing  vol¬ 
ume  available  in  a  porous  material;  the  Dilute  Spherical  Pores  (DSP) 
model  of  a  solid  with  a  low  density  of  spherical  pores  proposed  by 
Hasselman  [10]  and  referred  to  as  the  ‘non-linear’  semi-empirical 
model  by  Selcuk  and  Atkinson  [20];  and  the  theoretically  based 
Composite  Spheres  Model  (CSM)  proposed  by  Ramakrishnan  and 
Arunachalam  [17]  for  determining  the  bulk  elastic  modulus  of  a 
composite  material  composed  of  touching  spheres  of  one  or  more 
materials.  In  this  work  the  predictions  of  partially  porous  material 
properties  from  fully  dense  material  properties  has  been  shown 
to  agree  well  with  data  for  YSZ  and  YSZ  composites  [20].  Table  3 
presents  the  bulk  material  properties  of  the  composite  anode  as  a 
function  of  porosity  calculated  according  to  the  non-linear  (DSP) 
model  [20],  the  cathode  material  properties  are  from  Giraud  and 
Canel  [8]. 

Characteristic  strength  values  at  room  temperature  (25  °C)  and 
800  °C  are  taken  from  Selcuk  and  Atkinson  [20].  The  characteris¬ 
tic  strength  of  the  Ni/YSZ  anode  at  room  temperature  and  800  °C 
is  187  MPa  and  124  MPa  respectively,  for  the  YSZ  electrolyte  it  is 


Table  3 

Anode  and  cathode  bulk  material  properties  at  room  temperature  and  800  °C  as  a  function  of  porosity. 


Porosity  (%) 

25  °C 

800  °C 

Anode  (NiO-YSZ) 

Cathode  (LSM) 

Anode  (NiO-YSZ) 

Cathode  (LSM) 

E  (GPa) 

V 

E  (GPa) 

V 

E  (GPa) 

V 

£  (GPa) 

V 

0 

207.22 

0.328 

110 

0.36 

171.2 

0.331 

116.4 

0.36 

10 

161.5 

0.313 

- 

- 

133.4 

0.316 

- 

- 

20 

126.5 

0.301 

- 

- 

104.5 

0.304 

- 

- 

30 

99.0 

0.292 

41 

0.28 

81.8 

0.295 

43.4 

0.28 

40 

76.75 

0.283 

- 

- 

63.4 

0.286 

- 

- 

50 

58.4 

0.281 

— 

— 

48.3 

0.277 

— 

— 

Table  4 

Visco-plastic  material  properties. 


R.  Clague  et  al.  /  Journal  of  Power  Sources  210  (2012)  224-232 


227 


Property 

Temp  (K) 

Electrolyte  (YSZ) 

Cathode  (LSM) 

Anode  (Ni/YSZ) 

nc 

1273 

1.0  [4] 

1 .3  [  1 8  ] 

1.0  [9] 

Qc  (kj  mol-1) 

1273 

390-460 [4] 

405-521 [18] 

390-460  [9] 

A 

- 

40a 

5700b 

40c 

a  Calculated  from  Chokshi  [4]. 
b  Calculated  from  Routbort  et  al.  [18]. 
c  Calculated  from  Gutierrez  et  al.  [9]. 


232  MPa  and  154  MPa,  and  for  the  LSM  cathode  it  is  52  MPa  and 
75  MPa.  The  cathode  values  are  a  little  surprising  at  first  glance 
but  explained  by  the  material  becoming  more  ductile  at  800  °C  so 
more  plastic  deformation  is  tolerated  around  crack  tips,  leading  to 
a  higher  bulk  strength. 

2.2.  Visco-plasdc  material  behaviour 

Creep  is  the  term  used  to  describe  the  tendency  of  a  material 
to  relieve  internal  stress  over  time  by  deformation  at  stress  lev¬ 
els  below  the  material  yield  stress.  Creep  normally  occurs  when  a 
material  is  subjected  to  high  temperatures,  and  is  more  rapid  under 
high  stress.  The  stress  relieving  deformation  can  be  plastic  and  per¬ 
manent,  which  is  usually  the  case  with  metals  and  ceramics,  or 
elastic,  as  is  the  case  with  some  polymers  [6].  These  mechanisms  are 
defined  as  visco-plastic  and  visco-elastic  creep  respectively.  Metals 
typically  start  to  deform  visco-plastically  at  temperatures  around 
30-40%  of  their  melting  temperature,  with  ceramics  this  happens 
around  40-50%  of  their  melting  temperature  [1].  The  steady  state 
creep  behaviour  of  a  material  can  be  described  by  the  following 
relation  [1,6] 

e=/iexp(-|)  (i) 

where  6  is  the  steady-state  creep  rate,  A  is  a  material  constant,  Qc 
is  the  creep  activation  energy,  R  is  the  universal  gas  constant  and  T 
is  the  absolute  temperature.  The  material  constant  A,  has  complex 
dependencies  on  the  applied  stress  cr,  the  stress  exponent  nc,  the 
grain  size  L,  and  a  grain  size  exponent  ngs  and  has  the  form  given 
in  Eq.  (2) 


3.  Analytical  model  for  stress  analysis 

The  analytical  model  used  to  cross  check  the  predictions  of  the 
finite  element  model  is  based  on  original  work  by  Hsueh  [12]  where 
the  stress  distribution  in  the  nth  layer  of  a  multi-layer  is  calculated 
by  considering  the  other  n  -  1  layers  to  be  a  single,  composite  bar 
with  composite  material  properties  derived  by  averaging  the  mate¬ 
rial  properties  of  each  layer,  weighted  by  layer  thickness.  The  Hsueh 
model  allows  the  layers  to  have  individual  thermo-elastic  mate¬ 
rial  properties.  This  formulation  was  chosen  because  it  provides  an 
exact  solution  for  the  stress  in  the  layers  from  a  set  of  closed  form 
equations. 

A  multi-layer  strip  is  shown  schematically  in  Fig.  4  in  which  n 
layers  of  individual  thicknesses  q  are  laminated  together  such  that 
there  is  strain  continuity  at  layer  interfaces.  The  subscript  i  denotes 
the  layer  number  from  1  to  n,  the  material  properties  of  Young’s 
modulus  and  thermal  expansion  coefficient  are  denoted  by  E,-  and 
oi[  respectively. 

The  relationship  between  the  cumulative  thickness  h*  and  the 
individual  layer  thickness  q,  is  described  by 


hi  -  ytj  (i  =  1,  n)  (5) 

j=l 

The  Hsueh  approach  decomposes  the  total  strain  in  the  system  into 
a  uniform  component  of  an  effective  composite  bar,  and  a  bending 
component  induced  by  the  deviation  of  strain  in  a  layer  from  that 
in  the  composite  bar.  The  total  strain  in  the  system  is  therefore 
formulated  as  (6) 


A  =  crncL-nssp(  02)n°2  (2) 

Visco-plastic  material  behaviour  is  represented  in  Abaqus™  by  a 
power  law  creep  model  of  the  form  shown  in  Eq.  (3) 

£  =  Bonc  (3) 

where  6  is  the  strain  rate,  a  is  the  applied  stress,  nc  is  the  stress 
exponent  and  B  is  a  temperature  dependent  material  parameter 
normally  calculated  from  test  data  by  plotting  log  e  against  nc  log  a 
and  linearising  Eq.  (3)  by  taking  logs  of  both  sides  of  the  equation. 
Comparing  Eqs.  (1 )  and  (3)  and  assuming  no  dependence  on  partial 
pressure  of  oxygen  (i.e.  n02  =  0)  yields  the  form  of  B  as  Eq.  (4). 


6  =  c  +  — — —  (for  0  <  x3  <  hn)  (6) 

r 

where  c  is  the  uniform  strain  component  given  by  (7),  r  is  the  radius 
of  curvature  of  the  system  (8),  x3  is  the  through  thickness  position 
at  which  the  strain  is  calculated  and  tb  denotes  the  position  of  the 
bending  axis  (9)  in  direction  x3  taking  the  free  surface  of  layer  1 
as  x3  =  0.  The  bending  axis  is  a  definition  that  is  particular  to  the 
Hsueh  model  and  defined  as  the  line  in  the  cross  section  of  the 
system  where  the  bending  strain  component  is  zero.  This  is  not  the 
same  as  the  conventional  neutral  axis,  which  is  defined  as  the  axis 


B  =  L  exp  (-^)  (4) 

2.3.  Visco-plastic  material  properties 

Table  4  presents  the  visco-plastic  material  properties  used  in 
this  work.  As  the  material  properties  vary  significantly  with  grain 
size  it  is  critical  to  have  material  properties  that  are  derived  from 
appropriate  samples  with  representative  grain  sizes  and  material 
constants,  when  using  an  analysis  to  investigate  a  specific  problem. 
As  such  the  material  data  used  in  this  work  are  representative  only. 


Layer  n 

Layer i 

ti 

Layer  2 

^2 

Layer  1 

ti 

*3 


K-i 

h2 

h, 


0 


Fig.  4.  Definition  of  variables  in  analytical  multilayer  thermal  stress  model. 
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in  the  cross  section  of  the  system  where  there  are  no  longitudinal 
stresses  or  strains. 


AT  'y^EjtjOii 


c  = 


i=i 


i=l 


E\t^(2t\  +  3 tjj)  +  ^  +  6hi_\ti  +  2 1?  —  3tb(2h,_i  +  t,)] 

3[£i(c  -  ai  Ar)q2  -  ~  aiATX2hi- 1  +  fi)] 

n 

—E\t2  +  ^  ^£jf,(2ht_1  +  tj) 


f b  — 


i=2 


(7) 


(8) 


(9) 


layers,  contains  the  bending  axis,  which  changes  the  stress  distri¬ 
bution  in  the  anode  layer  so  that  part  of  it  now  experiences  tensile 
stress.  This  is  significant  because  fracture  of  the  ceramic  layers  will 
only  occur  under  tensile  stress,  thus  the  probability  of  failure  of  the 
structural  support  layer  containing  the  bending  axis  would  be  zero 
for  the  ‘planar’  case  in  Fig.  5  but  some  finite  value  for  the  ‘free’  case. 

4.  Finite  element  model  for  stress  analysis 

A  finite  element  model  was  generated  following  the  geometry 
shown  in  Fig.  1.  Using  the  C3D20  second  order  hexahedral  element 
from  the  Abaqus  element  library,  the  solution  was  demonstrated  to 
be  mesh  independent  when  the  number  of  second  order  elements 
was  greater  than  31616,  which  took  688  CPU  seconds  to  solve  a 
linear  static  stress  analysis.  All  analyses  used  a  computer  with  a 
Pentium  4,  2.8  GHz  CPU  with  1 .5  Gb  RAM. 

4.2.  Finite  element  analysis  validation 


The  in-plane  stresses  in  the  layers,  crf  are  related  to  the  in-plane 
strains  by  Eq.  (10) 

cri=Ei(€-asAT)  (10) 

If  a  plate  rather  than  a  strip  is  being  considered,  as  is  the  case 
with  this  SOFC  model  the  biaxial  Young’s  modulus  £'  should  be 
substituted  for  £*. 


where  v  is  Poisson’s  ratio  and  the  subscript  i  denotes  the  layer. 

The  mechanical  boundary  conditions  that  are  applied  to  a  stress 
model,  analytical  or  otherwise,  often  have  a  large  influence  on  the 
predictions  of  the  model.  In  this  analysis  the  SOFC  ceramic  mem¬ 
brane  has  been  simplified  to  a  three  layer  system  by  assuming  that 
the  electrode  active  layers  have  the  same  thermo-mechanical  prop¬ 
erties  as  their  counterpart  electrode  porous  layer.  In  order  to  asses 
two  extreme  sets  of  boundary  conditions  that  a  plate  might  expe¬ 
rience,  the  stresses  were  calculated  for  the  three  layer  laminated 
plate  when  it  is  constrained  to  remain  flat  (r  =  oo)  and  when  it  is  free 
to  bend  under  the  application  of  a  uniform  temperature  increase 
of  1000°C.  The  porous  anode  support  is  taken  as  layer  1,  the  elec¬ 
trolyte  as  layer  2  and  the  cathode  as  layer  3.  E[  and  are  taken  to  be 
uniform  and  constant  throughout  a  layer  according  to  Table  3.  Fig.  5 
presents  the  in-plane  stress  through  each  layer  of  the  simplified 
three  layer  SOFC  model. 

If  the  cell  is  constrained  to  remain  flat  the  stress  in  the  structural 
anode  layer  is  uniform  and  compressive.  When  the  cell  is  allowed 
to  bend  the  anode  layer,  being  substantially  thicker  than  the  other 


The  finite  element  model  was  cross  checked  by  considering  the 
case  of  a  three  layer  SOFC  heated  from  room  temperature  to  800  °C. 
The  stresses  predicted  by  the  finite  element  model  were  compared 
to  the  analytical  solution  for  two  limiting  cases;  the  first  ‘planar’ 
case  allowed  no  displacement  in  the  x3  direction  of  the  nodes  on 
the  x\  \x2  plane  at  x3  =0,  thereby  constraining  the  plate  to  remain 
flat  while  avoiding  over-constraint  by  allowing  in  plane  thermal 
expansion  and  through  thickness  expansion  in  x3 ;  the  second,  ‘free’ 
case  allowed  displacement  in  all  degrees  of  freedom,  in  effect  it 
is  an  unconstrained  plate.  Comparing  the  finite  element  analysis 
results  to  the  analytical  solution  (Fig.  6)  through  the  plate  thickness 
and  at  its  centre,  for  a  temperature  increase  of  800  °C  and  a  zero 
stress  condition  at  room  temperature,  shows  an  average  difference 
in  stress  between  the  numerical  and  analytical  predictions  of  4% 
for  both  cases.  The  difference  is  likely  to  be  due  to  edge  effects  in 
the  numerical  model  providing  stress  relief,  whereas  the  analytical 
model  assumes  the  domain  to  be  infinite  in  the  in-plane  dimension, 
thus  the  analytical  and  finite  models  are  considered  to  show  good 
agreement. 

4.2.  Residual  stress  model 

The  residual  stress  from  manufacture  was  modelled  by  select¬ 
ing  an  appropriate  temperature  at  which  the  SOFC  is  in  a  ‘zero 
stress  state’  on  cooling  from  sintering  and  using  this  as  the  ini¬ 
tial  conditions  for  the  analyses.  Fig.  7  shows  the  in-plane  principal 
stress  predicted  by  the  analytical  model  through  the  thickness 
of  the  ceramic  membrane  for  two  different  cases:  residual  stress 
neglected;  residual  stress  included.  The  analysis  assumes  a  zero 
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Fig.  7.  Comparison  of  stress  predictions  including  residual  stress  and  without  resid-  Fig.  8.  Boundary  condition  comparison, 

ual  stress. 


stress  temperature  of  25  °C  for  the  first  case  and  1010°C  for  the 
second  following  Atkinson  and  Selcuk  [3],  both  have  a  uniform 
operating  temperature  of  800  °C  imposed. 

Fig.  7  shows  clearly  that  the  stress  state  of  the  ceramic  mem¬ 
brane  at  800  °C  is  affected  a  great  deal  by  the  inclusion  of  residual 
stress,  which  changes  the  sign  and  magnitude  of  the  predicted  in¬ 
plane  principal  stresses  in  the  cathode  and  electrolyte  from  mean 
values  of  20  MPa  and  308  MPa  when  residual  stress  is  neglected 
to  -5.3  MPa  and  -83  MPa  respectively  when  residual  stress  is 
included.  The  stress  levels  with  no  residual  stress  are  very  signif¬ 
icant  compared  to  the  characteristic  strengths  of  the  cathode  and 
electrolyte  of  75  MPa  and  154  MPa  (at  800  °C),  the  probability  of 
failure  is  likely  to  be  1  or  nearly  1 .  When  residual  stress  is  included 
the  stresses  are  compressive  and  therefore  unlikely  to  cause  brittle 
fracture  of  the  material. 

The  in-plane  principal  stress  in  the  anode  support  layer  also 
changes  sign  and  magnitude  when  residual  stresses  are  included. 
At  x3  =  0  it  changes  from  9.1  MPa  without  residual  stress  to 
-2.66  MPa  when  residual  stress  is  included,  and  at  x3  =  610p>m 
(the  electrolyte-anode  interface)  it  changes  from  -17.3  MPa  with¬ 
out  residual  stress  included  to  5.04  MPa  with.  Although  the  tensile 
stress  levels  in  the  anode  are  small  compared  to  the  anode  char¬ 
acteristic  stress  of  124  MPa  (at  800  °C)  there  will  still  be  a  finite 
probability  of  failure. 

4.3.  Mechanical  boundary  conditions 

Three  different  sets  of  mechanical  boundary  conditions  were 
defined  for  the  finite  element  model,  these  are;  free  to  bend;  con¬ 
strained  to  remain  flat;  bonded  around  its  edge  (represented  by 
surface  Si  in  Fig.  1 )  into  a  test  housing  at  room  temperature  after 
cooling  from  sintering.  These  are  summarised  in  Table  5.  A  com¬ 
parison  of  the  peak  predicted  stress  in  each  layer  of  the  fuel  cell 
for  each  of  the  three  boundary  conditions  is  shown  in  Fig.  8.  It  can 
be  seen  that  the  different  mechanical  boundary  conditions  have  a 


Table  5 

Definition  of  boundary  condition  sets. 


Name 

Definition 

Description 

BC1 

At  X3  =  0,  no 
displacement  in  X3 
direction 

Membrane  remains  planar 

BC2 

No  restraint 

Membrane  free  to  bend 

BC3 

Surface  SI, 

X3  =  fixed  after 
cooling  from 

1010°C 

Representative  of  test  cell 

significant  impact  on  the  calculated  stress  distribution.  Boundary 
condition  set  BC3  produces  the  highest  stress  in  the  anode  layer, 
which  is  located  in  the  area  where  the  anode  is  bonded  to  the  hous¬ 
ing.  Set  BC2,  where  the  SOFC  is  constrained  to  remain  flat,  gives  the 
highest  stress  in  the  electrolyte  and  cathode  layers  as  there  is  no 
strain  relief  due  to  the  curvature  of  the  SOFC. 

4.4.  Redox  model 

The  change  in  anode  properties  on  chemical  reduction  is  mod¬ 
elled  by  defining  material  properties  that  are  dependent  on  a  field 
variable  representing  the  porosity,  as  given  in  Table  3.  The  poros¬ 
ity  field  variable  is  defined  using  the  Abaqus  *FIELD  keyword  and 
allowed  to  evolve  linearly  over  one  analysis  step  from  an  initial 
value  of  20%,  representing  NiO/YSZ,  to  a  final  value  of  40%,  rep¬ 
resenting  Ni/YSZ.  The  evolution  of  the  maximum  and  minimum 
in-plane  stresses  as  the  anode  is  chemically  reduced  is  shown  in 
Fig.  9. 

The  change  in  the  layer  in-plane  principal  layer  stresses  pre¬ 
dicted  by  the  finite  element  model  are  compared  with  the 
predictions  of  the  analytical  model  in  Table  6  in  order  to  validate 
the  method  used  in  the  finite  element  model.  The  results  are  simi¬ 
lar  for  both  the  analytical  and  finite  element  models,  and  show  that 
a  reduction  in  peak  stress  of  more  than  15%  is  experienced  in  the 
cathode  when  the  anode  undergoes  chemical  reduction.  The  dif¬ 
ference  between  the  finite  element  analysis  and  analytical  model 
is  once  again  ascribed  to  edge  effects  in  the  finite  element  model,  a 


Normalised  Analysis  Step  Time 


Fig.  9.  Stress  evolution  in  layers  on  reduction  of  anode. 
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Table  6 

Comparison  of  stress  predictions  on  anode  reduction,  FEA  and  analytical  models 
(MPa). 


FEA  (MPa) 

Analytical  (MPa) 

Ox 

Red 

%  Diff 

Ox 

Red 

%  Diff 

Anode  (top) 

5.0 

4.6 

-8.7 

5.2 

4.8 

-7.5 

Anode  (bttm) 

-2.6 

-2.3 

-8.9 

-2.7 

-2.5 

-7.5 

E’lyte 

-87.2 

-83.3 

-4.4 

-86.7 

-83.4 

-3.8 

Cathode 

-6.3 

-5.2 

-17.7 

-6.2 

-5.2 

-15.5 

consequence  of  the  finite  element  model  domain  being  finite  rather 
than  infinite  as  in  the  analytical  model. 

4.5.  Visco-plasticity  model 


Using  the  model  of  visco-plasticity  described  in  Section  2.2  the 
stress  evolution  over  time  was  analysed.  In  this  analysis  the  cell 
was  free  to  bend  and  the  zero  stress  condition  was  taken  as  25  °C. 
Following  the  discussion  in  Section  2.2  the  visco-plastic  properties 
of  the  Ni/YSZ  support  layer  were  considered  to  be  those  of  pure 
YSZ.  Fig.  10  shows  the  time  dependent  stress  evolution  predicted 
by  the  finite  element  model  at  a  point  in  the  centre  of  the  fuel  cell 
and  x3  =  0  (the  bottom  of  the  anode  layer),  at  800,  900  and  1 000  °C 
over  10,000  h. 

It  can  be  seen  that  the  operating  temperature  of  the  fuel  cell 
makes  a  large  difference  to  the  rate  at  which  stresses  are  relieved 
by  visco-plastic  creep.  At  1000°C  the  stress  is  completely  relieved 
after  15  x  106  s,  whereas  at  800  °C  the  stress  has  only  reduced  by 
approximately  2%.  This  analysis  suggests  that  when  operating  at 
1 000  °C  SOFCs  significant  stress  relief  by  visco-plastic  creep  is  likely 
to  occur  during  the  operating  lifetime  of  the  SOFC.  It  is  also  likely 
that  creep  under  gravity  will  occur  at  these  temperatures. 

4.6.  Temperature  boundary  conditions 

A  uniform  temperature  distribution  has  been  assumed  above 
to  allow  an  analytical  model  to  be  developed  and  used  to  vali¬ 
date,  to  some  degree,  a  finite  element  model.  A  more  sophisticated 
approach  is  to  use  a  CFD  model  to  provide  nodal  temperature 
distribution  for  use  in  the  FE  model.  The  finite  element  nodal  tem¬ 
peratures  were  interpolated  from  the  results  of  a  CFD  simulation.  In 
this  way  the  CFD  and  FE  models  are  coupled  and  the  need  to  make 
simplifying  assumptions  like  constant  through  thickness  temper¬ 
ature  is  eliminated.  Fig.  1 1  shows  the  temperature  distribution  at 


Fig.  11.  Temperature  distribution  at  an  electrical  load  of  300  mA cm-2. 


Table  7 

Definition  of  temperature  loadcases. 


Name 

Definition 

Description 

O 

O 

kT 

Uniform  1010°C 

Sintering  temperature 

T25 

Uniform  25  °C 

Room  temperature 

^800 

Uniform  800 °C 

Nominal  furnace  temperature 

Tcfd3 

From  CFD 

Temp  distribution  at  300  mAcirr2 

an  average  electrical  load  of  300  mA cm-2.  The  definitions  of  the 
temperature  loadcases  referred  to  later  on  are  shown  in  Table  7. 

The  results  of  the  finite  element  analysis  of  the  maximum  prin¬ 
cipal  stress  distribution  throughout  the  SOFC  ceramic  membrane 
for  the  temperature  boundary  conditions  T80o  and  Tcfd3  are  shown 
in  Figs.  12  and  13  respectively.  In  these  figures  only  half  of  the  fuel 
cell  is  shown,  the  left  hand  image  is  a  view  on  the  top  of  the  anode 
layer,  the  central  image  is  a  top  view  on  the  electrolyte  layer,  and 
the  right  hand  image  is  a  top  view  on  the  cathode  layer.  Below  each 
of  the  top  views  is  a  perpendicular  section  through  the  layer.  In 
these  analyses  the  zero  stress  temperature  is  taken  to  be  1010°C. 

The  largest  variation  in  maximum  principal  stress  through  a 
layer  thickness  is  seen  in  the  anode,  primarily  because  it  is  the  thick¬ 
est  layer  and  contains  the  neutral  axis  of  the  system.  For  the  CFD 
temperature  distributions  the  thermal  stress  distributions  caused 
by  the  thermal  gradients  induced  by  the  channels  and  plena  of  the 
fuel  cell  housing  is  clearly  visible  by  comparing  Fig.  12  with  Fig.  13. 
The  thermal  gradients  due  to  the  electrochemical  reactions  in  the 


Fig.  10.  Visco-plastic  stress  relaxation  over  1 0,000  h. 


Section  X-X.  Anode  Section  X-X,  Electrolyte  Section  X-X,  Cathode 


Fig.  12.  Maximum  principal  stress  distribution  in  SOFC  at  uniform  temperature  of 
800  °C. 
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Fig.  13.  Maximum  principal  stress  distribution  in  S0FC  at  an  electrical  load  of 

300mAcm-2.  Time  (hrs) 


Fig.  15.  Maximum  principal  stress  evolution  in  layers  over  end  of  duty  cycle  using 
fuel  cell  induce  significant  additional  stresses  when  compared  to  CFD  derived  boundary  conditions. 

the  uniform  temperature  case. 


5.  Stress  evolution  over  a  complete  duty  cycle 

The  analyses  so  far  have  considered  separate  mechanisms  of 
stress  generation  and  relaxation  that  are  active  in  the  SOFC  from 
manufacture  to  the  start  of  electricity  generation.  To  determine  a 
stress  history  for  the  fuel  cell  ceramic  membrane,  all  of  these  mod¬ 
els  must  be  combined  into  one  finite  element  analysis  in  the  time 
domain,  with  appropriate  mechanical  restraints  applied. 

5.2.  Stress  evolution  over  duty  cycle 

The  in-plane  stress  evolution  in  the  anode  layer  over  time  is 
shown  in  Fig.  14  for  the  SOFC  duty  cycle  defined  in  Fig.  3  where 
the  initial,  zero  stress  temperature  is  taken  as  1 01 0  °C.  As  the  SOFC 
cools  from  sintering  between  time  T=  0  and  T=  10  h  the  peak  stress 
in  each  layer  increases.  At  room  temperature  the  cell  is  bonded 
into  the  housing,  modelled  by  restraining  the  X3  degree  of  free¬ 
dom  of  the  nodes  on  surface  Si.  As  the  cell  is  heated  up  to  800 °C 
from  T=  10  and  T=  1 7  h,  the  stress  levels  drop  as  the  membrane  is 
now  closer  to  its  sintering  temperature.  Between  T=  1 7  and  T=  20  h 
the  anode  material  is  chemically  reduced  from  NiO  to  Ni  leading 


to  a  drop  in  the  bulk  elastic  modulus  of  the  anode.  This  leads  to  a 
change  in  the  stress  distribution  such  that  the  peak  stress  in  the 
anode  and  electrolyte  layers  drops,  while  that  in  the  cathode  layer 
increases.  The  next  step  in  the  duty  cycle,  from  T=  20  to  T=  2000  h, 
is  to  draw  current  and  as  a  temperature  distribution  representing  a 
low  current  density  of  70  mA  cm-2  is  applied  the  peak  stress  in  the 
anode  and  cathode  layers  increases.  During  this  time  the  stresses 
are  gradually  relieved  by  visco-plastic  creep  but  as  the  operating 
temperature  is  approximately  800  °C  the  effect  is  small.  At  T=  2000 
to  T=4000  h  a  temperature  distribution  due  to  an  electrical  load  of 
300  mA  cm-2  is  applied,  leading  to  a  sharp  increase  in  peak  stress  in 
all  layers  and  from  T = 4000  to  T=  6000  h  a  temperature  distribution 
due  to  a  mid-range  electrical  load  of  200  mAcnrr2  is  applied,  with 
a  corresponding  drop  in  peak  stress  levels.  The  final  stress  state  of 
the  layers  can  be  seen  in  detail  in  Fig.  15,  which  shows  the  stress 
evolution  over  the  final  part  of  the  duty  cycle,  i.e.  cooling  to  room 
temperature.  The  final  stress  state  of  the  SOFC  layers  at  room  tem¬ 
perature  after  being  subjected  to  the  duty  cycle  is  very  different 
to  that  at  room  temperature  when  T=  lOh,  for  example  the  peak 
stress  in  the  electrolyte  atT=10hisll  MPa,  whereas  at  6030 h  it  is 
17  MPa. 


Time  (hrs) 


Fig.  14.  Maximum  principal  stress  evolution  in  SOFC  layers  over  duty  cycle. 


6.  Conclusions  and  closing  remarks 

In  this  paper  models  of  mechanisms  for  generating  and  relieving 
stress  in  the  ceramic  membrane  of  SOFC  have  been  developed: 

•  thermo-mechanical  stress  due  to  differential  thermal  expansion 
or  temperature  gradients; 

•  residual  stress; 

•  stress  due  to  membrane  mechanical  restraint; 

•  stress  change  due  to  anode  reduction; 

•  stress  relaxation  by  visco-plastic  creep. 

The  inclusion  of  residual  stress  has  a  significant  effect  on  the 
stress  magnitude  and  distribution  in  the  SOFC  at  operating  condi¬ 
tions,  and  is  considered  the  most  important  aspect  to  include.  When 
cooling  from  sintering  to  room  temperature  and  heating  from  room 
temperature  to  operating  temperature,  stresses  are  generated  due 
to  differential  thermal  expansion  of  the  layers  that  are  similar  to 
the  characteristic  stress,  indicating  that  there  is  a  high  probability 
of  failure  during  these  phases  of  the  life  cycle.  The  redistribution  of 
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stress  caused  by  the  change  in  material  properties  when  the  anode 
is  chemically  reduced  from  NiO/YSZ  to  Ni/YSZ  is  of  smaller  magni¬ 
tude  than  the  inclusion  of  residual  stress  but  still  significant  with 
the  largest  changes  in  predicted  stress  being  of  the  order  of  15%. 
The  inclusion  of  a  model  for  visco-plastic  creep  of  YSZ  has  indi¬ 
cated  that  significant  stress  relief  will  occur  over  SOFC  operating 
timescales.  This  has  implications  for  sealing  methods  that  rely  on 
compression  of  the  cells  as  the  area  under  the  seal  will  tend  to 
relax  over  time.  It  also  implies  that  any  SOFC  with  a  support  struc¬ 
ture  based  on  YSZ  that  is  loaded,  even  by  gravity,  will  deform  over 
timescales  of  a  few  thousand  hours  at  SOFC  operating  tempera¬ 
tures.  The  difference  in  stress  distribution  resulting  from  the  T800 
and  Tcfd3  temperature  distributions  has  been  shown  to  be  impor¬ 
tant  with  all  three  layers  showing  significant  tensile  stresses  for  the 
CFD  derived  temperature  distribution,  while  for  the  uniform  tem¬ 
perature  approximation  the  stresses  in  the  cathode  and  electrolyte 
layers  are  predominantly  compressive  and  therefore  unlikely  to 
cause  failure.  Comparing  the  a  1 1  (=  o2  2 )  stress  of  the  uniform  800  °C 
case  with  the  maximum  principal  stress  of  the  CFD  derived  Tcfd3 
boundary  conditions  it  can  be  seen  that  the  tensile  stress  is  much 
higher  for  the  latter  condition,  indicating  that  the  temperature 
distribution  generated  by  the  electrochemical  reactions  generate 
significant  additional  stresses. 
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